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ABSTRACT.  A  sensitivity  methodology  for  nonlinear  delay  systems  arising  in  one  class  of  cellular  HIV  infection  models  is  pre¬ 
sented.  Theoretical  foundations  for  a  typical  sensitivity  investigation  and  illustrative  computations  are  given. 


1.  Introduction 

1.1.  Background.  Over  the  past  several  years,  the  use  of  mathematical  models  as  an  aid  in  understanding  features  of 
HIV  and  other  virus  infection  dynamics  has  been  substantial.  Several  papers  published  in  the  mid  nineties  provided  strong 
evidence  for  the  high  rate  of  HIV-1  replication  and  clearance  in  infected  individuals  [16,  37,  51].  By  the  end  of  the  decade, 
the  general  consensus  was  that  in  vivo,  on  the  order  of  10 10  virions  are  assembled  and  cleared  every  day  [25,  35,  39], 
In  many  of  these  papers,  the  viral  clearance  rate  c  was  identified  by  modeling  the  disease  pathogenesis  with  a  system  of 
deterministic  differential  equations,  numerically  calculating  a  solution,  and  then  fitting  the  results  with  experimental  data 
(using  a  nonlinear  least  squares  (NLS)  approach),  e.g.,  see  [35,  37,  39].  The  existence  of  such  a  high  replication/clearance 
rate  implies  a  high  mutation  rate,  thus  indicating  that  pharmacological  mono-therapy  will  ultimately  fail,  since  the  virus  can 
rapidly  manifest  a  resistance  to  any  one  drug.  More  importantly,  this  knowledge  directly  contributed  to  the  current  practice 
of  simultaneously  administering  multiple  drugs  to  HIV  positive  individuals  in  an  effort  to  counteract  the  high  mutation  rate 
of  the  virus. 

Following  its  success  in  helping  to  identify  this  significant  feature  of  the  HIV  pathogenesis,  the  use  of  mathematical 
modeling  and  parameter  identification  in  the  study  of  HIV  experienced  a  dramatic  increase.  In  particular,  in  the  wake  of 
the  publication  of  [37],  there  were  papers  covering  everything  from  additional  and/or  alternative  compartment  formulations 
[7,  21,  27,  28,  32,  38,  41,  52,  56,  57]  to  arguments  for  and  against  the  use  of  delay  differential  equations  in  modeling  the 
eclipse  phase  [12,  13,  15,  24,  26,  29,  30,  31]  (including  those  that  addressed  the  solution  stability  [30,  31,  44]).  Moreover, 
in  the  context  of  delay  equations,  many  of  these  papers  focused  heavily  on  the  inter-relationship  between  the  parameters 
describing  the  drug  efficacy  //,  the  length  of  the  eclipse  phase  r,  the  infected  T-cell  death  rate  S,  and  the  virion  clearance  rate 
c  [12,  15,  24,  26,  29,  30,  31,  44].  The  purpose  of  this  paper  is  to  illustrate  our  approach,  which  allows  one  to  develop  new 
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insights  into  HIV  pathogenesis  by  utilizing  a  mathematical  tool  not  typically  associated  with  conventional  NLS  techniques. 
Indeed,  there  is  a  precedence  for  this  approach,  as  is  evidenced  by  previous  papers  within  the  HIV  modeling  literature 
that  make  use  of  stochastic  analysis  and  inference  [18,  45,  46,  49,  54,  58],  control  theory  [19,  53],  and  nonlinear  analysis 
[14,  50],  Note  that  the  above  survey  is  not  intended  to  be  comprehensive,  as  there  already  exist  thorough  reviews  of  the  field 
presented  in  [33,  34,  36], 

For  any  system  of  differential  equations  designed  to  model  real  world  phenomena,  whether  it  be  biological,  chemical,  or 
physical,  a  common  goal  is  to  understand  the  manner  in  which  the  system’s  constitutive  parameters  influence  its  solution. 
These  parameters  (such  as  S  above)  are  designed  to  correspond  to  aspects  of  the  phenomena  under  investigation  (e.g.,  HIV 
pathogenesis),  and  thus  it  is  desirable  to  predict  how  changes  in  the  parameters  will  affect  the  solution.  Indeed,  there  are 
several  papers  in  the  HIV  modeling  field  which  focus  heavily  on  the  topic  (good  examples  include  [41,  43]).  One  way  to 
address  this  question  is  to  perform  a  sensitivity  analysis,  a  mathematical  tool  developed  in  the  context  of  modern  control 
theory  and  commonly  used  in  mechanical,  aerospace,  electrical,  and  structural  engineering.  The  precursor  of  this  technique 
can  be  traced  back  to  an  1833  electrostatics  experiment  designed  to  measure  the  inductance  of  certain  metals  [8].  However, 
significant  activity  in  this  area  only  arose  in  the  middle  part  of  this  century,  concomitant  with  the  development  of  modern 
control  theory  in  the  late  1930’s.  In  our  analysis,  we  will  employ  the  semirelative  sensitivity  function,  though  there  are  other 
possibilities,  such  as  the  logarithmic  sensitivity  function  advocated  by  Bode  in  his  book  on  electrical  network  analysis  [5], 
We  direct  the  interested  reader  to  the  following  introductory  texts  [10,  1 1],  advanced  texts  [20,  47,  48,  55],  and  surveys  of 
the  field  [1,9].  We  also  note  that  the  sensitivity  analysis  described  in  this  paper  should  not  be  confused  with  the  statistical 
technique  of  the  same  name  and  based  on  Latin  Hypercube  Sampling  [4,  17]. 


1 .2.  Approach.  The  first  step  in  the  sensitivity  analysis  is  to  derive  the  sensitivity  equations  by  formally  taking  derivatives 
(with  respect  to  a  parameter  of  interest)  on  both  sides  of  the  original  equation(s).  The  solution  to  this  new  system  (assuming 
for  the  moment  it  is  well-posed)  contains  information  regarding  the  sensitivity  of  the  original  system  to  perturbations  in 
the  chosen  parameter  (around  some  a  priori  fixed  value  of  that  parameter).  Hereafter  we  will  refer  to  the  solution  to  the 
sensitivity  equations  as  a  sensitivity  function. 

To  illustrate  the  sensitivity  procedure,  we  will  examine  an  HIV  population  system  with  compartments  described  in  [3], 
summarized  in  Table  1,  and  denoted  by  the  vector  x  =  ( V,A,C,T)t .  In  this  case  (see  [3]  for  details),  our  system  of 
distributed  delay  differential  equations  is 


(1.1) 


x(t)  =  L(x(t),xt)  +  fi(x(t))  +  /2(f)  for 0  <  t  <  tf 
(;c(0),  Xq)  =  ($(0),$)  G  M4  x  0;  R4)  , 
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where  r  and  tf  are  finite,  x t  denotes  the  function  6  ha  x{t  +  #),  0  £  [— r,  0],  and  for  (77,  <p)  £  l4  x  ^{—r,  0;  IR4), 


L(V,4>)  = 


—c  0  no  0 

0  rv  —5a  0  0 

0  0  rv  —  5c  0 

0  0  0  r„  —  5„ 


il  +  nA[5(  1,2)]  {4j4j  /  <pU>)dP\  {()) 

J  —r 


+7([^(3,2)](4j4)  -  [<5(2,2)]  (4,4)^  [  <t)(9)dP-2(d) 

’  J  —r 


fiOi) 


hit) 


-mm 

-siY)i=27ii)m  +  mm 
-<KE  t=2  ndm 
_  di.yj  -  pmm 
[0, 0, 0,  s]T  ,  0  <  t  <  tf  . 


Here  the  compartments  in  x  and  the  parameters  (including  the  probability  distributions  P 1,  Pi)  given  in  the  vector  q  = 
(c,rv,ru,nA,nc,5,5A,5c,5Ui7iP,  Pi,P2,  S)T  are  all  described  in  [3].  A  full  and  thorough  sensitivity  analysis  could 
include  not  only  derivatives  with  respect  to  the  scalar  parameters  (e.g.,  7  or  5  a X  but  also  Frechet  derivatives  with  respect 
to  the  delay  distributions  (e.g.,  P\  or  Pi).  The  following  sections  include  discussions  regarding  the  well-posedness  of  the 
sensitivity  equations  and  an  example  numerical  simulation  as  well  as  an  interpretation  of  the  results. 


Notation 

Description 

V 

Infectious  viral  population 

A 

Acutely  infected  cells 

C 

Chronically  infected  cells 

T 

Uninfected  or  target  cells 

Table  1.  in  vitro  model  compartments 


2.  Theory 

For  those  interested  in  the  mathematical  considerations,  this  section  contains  theoretical  foundations  for  the  well- 
posedness  of  the  sensitivity  equations.  While  the  results  presented  here  are  important  because  they  legitimize  our  study 
of  these  equations,  understanding  the  techniques  in  the  proofs  are  not  essential  to  appreciating  the  simulations  and  results 
presented  in  Section  3.  Therefore,  those  readers  who  wish  to  skip  the  details  in  this  section  may  do  so  with  little  loss  in 
understanding  the  formal  aspects  of  sensitivity  analyses. 

For  our  illustrative  discussions,  we  will  only  consider  distributions  Pi ,  Pi  that  are  both  differentiable  and  parameteriz- 
able  by  a  mean  p  and  a  standard  deviation  a  (i.e.,  for  i  =  1, 2,  pi(0)  =  J 'gPi(9 )  and  Pj(0)  =  Pj(6 ;  p,:,  cr,:)  for#  G  [— r,  0]). 
Moreover,  we  further  assume  that  the  resulting  densities  p,  are  tf1  in  p.j  and  a-,,  respectively.  To  illustrate  a  sensitivity 
analysis,  let  us  fix  the  forms  of  the  distributions  Pi,  Pi  and  consider  for  t  £  [—r.  tf],  the  derivative  of  x{t;  p  1 )  with  respect 
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to  /i]  (where  p\  is  the  parameter  corresponding  to  the  mean  of  pi).  If  we  let  (rp  <p)  G  IR4  x  (to{—  r,  0;  IR4),  t  G  [0,f/], 
pi  >  0,  then  from  results  established  in  [6],  we  note  that  lf(t,  r],  4>,p\)  =  L(rj,  <j> ;  p i)  +  fiff)  +  /2(f)  is  <^’1  in  t,  77,  (p,  and 
pi  under  smoothness  assumptions  (detailed  in  [3])  on  JL  L,  /  1,  and  /o.  For  our  specific  case,  to  prove  that  the  derivative 
of  x  with  respect  to  p  1  exists  and  is  continuous  in  /,  we  will  make  use  of  the  following  lemma. 

Lemma  2.1.  There  exists  a  solution  to 

^  y(i)  =  pi(x(t;pi);p(f)) +p2(pi;t/i) +P3(.'Ci(pi),pi;  1)  forO<t<tf 

(1/(0),  j/o)  =  m^GM'xff-r.O;!4), 

forx(t\Hi)  the  solution  to  (1.1),  and  where  for  p,  £  G  IR  r\,  (  G  IR4,  <p,ip  G  <^7(— r,  0;  IR4), 

1° 

92{/J;ip)  =  nA[S( i,2)]{44)  /  ip{9)p1{9;p,a1)d9 

+t(  [<5(3,2)]  {4_4)  -  [<5(2,2)]  (4j4))  /  ip{9)p-2 {9;  pa ,  cr-2)d9 

’  J—  r 

S3(4p;£)  =  'k4  [<5(1,2)]  {44)  y  o(0)(-^-i>{(0:ii.(T{))(Z)<IO . 

and  where 


-c  -  p/p 

0 

nc 

-pm 

p/p  r 

V  -  Sa  -  <5 (2/p  +  773  +  ?p) 

— <5/72 

— <5/72  +  pm 

0 

—<5/73 

-  <5c  -  <5  (  772  +  2/73  +  zp) 

— <5/73 

-p/p 

—<5/74 

— <5/74  L,  - 

■  <5U  -  <5(r72  +  773  +  2/p)  -  77/p 

Proof:  On  the  right  side  of  (2.1),  the  function  pi  +  <72  +  </3  satisfies  both  the  differentiability  condition  (Lemma  4.1) 
and  the  global  Lipschitz  condition  (Lemma  4.2)  from  [3],  Following  the  reasoning  in  the  proof  of  Theorem  4.5  in  the  same 
reference,  by  defining  a  convergent  sequence  of  successive  approximations,  it  can  then  easily  be  shown  that  a  solution  exists 
and  is  unique. 

Remark  2.2.  Note  that  Lemma  2. 1  guarantees  the  existence  of  a  solution  to  a  system  of  equations  with  a  general  initial 
condition  Tc  Recall  that  in  equation  (1.1),  the  initial  condition  <I>  is  independent  of  p  1  and  thus  the  next  step  will  be  to  argue 
that  system  (2.1)  combined  with  the  trivial  initial  condition  T<  =  0  comprises  the  sensitivity  equations. 

Theorem  2.3.  For  the  solution  x  of  (1.1),  x  has  a  derivative  with  respect  to  the  parameter  p  1  and  for  p  \  =  p  >  0,  this 
derivative  v(t.)  =  p)  satisfies  (2.1 )  with  the  initial  condition  ( (0) ,  ik)  =  (0,  0)  G  M4  x  ^(—r,  0;  IR4). 

Proof:  To  prove  the  existence  of  a  derivative  of  x  with  respect  to  the  parameter  p  1,  we  fix  p  \  at  p  >  0,  let  e  G  IR  be  a 
perturbation  of  p,  and  for  all  t  G  [— r,  t /],  define 


h(t,  p,  e)  =  x(t\ p  +  e)  —  x(t ;  p) . 
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The  overall  structure  of  the  proof  is  thus  to  show  that 


=  lim  LLA 

on i  H->o  e 


exists  and  is  continuous  for  t  G  [— r,  t /].  We  begin  by  considering  h 

h(t,n,e)  =  /  {^r{s,x{s-,n  +  s),xs(n  +  e),n  +  s)  -  ^{s,x{s\n),xs{n)^)}ds 

Jo 

=  /  n  +  +  c) ,  n  +  e)  —  J?(s,x(s;  n),xs(n  +  e),  ft  +  £) 

Jo 


+&{s,x{s]  n) , xs in  +  e),n  +  e)  -  &  ( s,x(s ;  p.) ,  (//) ,  n  +  e) 

+.?(.s.t(.s; //)../• .,(//).//  +  •  )  ,'?(«. x (, s: //)../•..(//).// ) }  (/.s  . 

According  to  the  Mean  Value  Theorem  [23],  we  have 

h(t,fj,e )  =  /  /  {Dx.^(s,x(s-,ij)  +  s'h(s,n,£),xs(n  +  s),n  +  e)(h{s,n,£)) 

Jo  Jo 

+DXt^{s,x{s;  /j),xs{d)  +  *'//„(//, s).//  +  y)(hJn,e)) 

+DIXl^'(s,  x(s ;  n)iXs(n),  n  +  s,£)(£)}  ds'ds , 

where  DXr,'P,  Dlri A  are  the  Frechet  derivatives  of  with  respect  to  its  second,  third,  and  fourth  arguments 

respectively.  Since  is  V' 1  in  all  its  arguments,  we  then  know  that  there  exists  linear  functions  A  s,  A~,  ,  A®,  _  (each 
parameterized  by  s '  and  e)  such  that 

h{t,ii,e)  =  {Dx^(s,x{s\fi),x:s{fi),n)(h(s,fj„£))  +  ^l',e(h{s,n,£)) 

Jo  Jo 

+  Dxr?(«.,x(»:  n),xs{p),n)(hs{n,  e))  +  A2s,  s(hs(n,  e)) 

+Dfll.^(s,x(s;n),xs(n),n)(£)  +  A^e)}  ds'ds  , 

where  |A4,  J,  |A^,  J,  |A^,  s\  — >  0  uniformly  in  .s'  as  >  — >  0.  Thus  for  s  G  [0,  f],  v,  £  G  US,  i],  (  G  M4,  <p,^  G  ^(—r,  0;  M4), 
and  <?i ,  <72 ,  <)?,  as  defined  in  Lemma  2.1, 


Then  the  equation  for  h  is 


Dx^(s,f],(j),v)(0 
DXt&{s,r],<j),v){ip) 
!),,  ./(s.i/.o.  i/)(0 


9i(m0 
HA''-  <■") 
9o(<P,v]0  ■ 
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Jo 

+  /  /  +  A2s,t£(hs(n,e))  +  A^(e)}  ds'ds. 

Jo  Jo 

Moreover,  since  gi,  g3,  g3  are  all  linear  in  their  last  arguments,  the  equation  for  h  can  be  used  to  obtain 


<  max 

rG[i-r,t]  . 


+ 


J  |  gi(x(s-,  fj,);  •)  +  J  \A],£\ds' 


\h(s,fi,£)\ 


g-An-  1+  /  \^l,e\ds' 

Jo 

+  93(Xs(f*),fJ,\  ■)  +  j  \Al,s\ds'  |  : 

where  1 1  •  |  j  is  the  oo-norm  on  [t  —  r,  t].  Thus  for  constants  K i ,  Kn  >  0,  we  know  that 

\\ht(fJi,e)\\  <  K\  f  \\hs(fi,s}\\ds  + K2tf\e\i 

Jo 

and  a  simple  application  of  Gronwall’s  inequality  implies  that 


(2.2) 


IIM^.e)ll  <  K-2  \tf\ |c|exp(A'1f/). 


which  will  be  useful  in  the  next  step. 

Now,  if  we  divide  both  sides  of  the  equation  for  h  by  e  so  that 

h(t,fi,£ ) 


gi{x(s,g)\  +  g-2(g-  +  g3{xs(g),  g:,  -)  j  ds 

£  £  £  ) 

t 


+ 


Jo  j  {A^e(^M)  +  A^(^l£l)  +  A^(£)}  ds'ds. 


we  note  that  the  form  of  the  integrand  is  strikingly  similar  to  the  right  side  of  the  equation  in  (2.1).  For  equation  (2.1), 
we  denote  the  solution  generated  using  g3  =  g  and  initial  condition  ('I'(O),  <&)  =  (0,  0)  G  M4  x  <^’(— r,  0;  R4)  as  v(t)  for 
t  G  [— r,  tf].  Moreover,  we  claim  that  this  solution  v  is  equal  to  the  limit  of  /i/e  as  Je  — >•  0. 

By  Lemma  2.1,  we  know  that  v  exists  and  is  continuous  for  t  G  [— r,  t /].  Clearly  v  and  h/%  are  identically  zero  for 
t  G  [— r,  0]  for  all  e  >  0,  and  thus  we  consider  for  t  G  [0,  t /] 


j{t)  - 


h(t,  g,  e) 


< 


Vt 


ht(g,£) 


<  max 

r£.[t—r,t\ 


n,  ,  x  ,  x  h(s,g,£ )  hs(g,£) 

g1(x(s,g);v(s) - 7 - )  +  g-2(g-,vs - 7 - J 


£ 

T  pi 


+g3{xs{g),g\l-^)^ds  -  j  J  j  A], 


+A2^(MM)+A3(j£)|  dTds 
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<  max 
re[t— r,t] 


n 

0  l 


\gi(x(s,fi);  -)| 


;(s)  - 


h(s,fi,,£) 


+  •  )l 


V*  - 


hs(fJ,,e) 


+ 


j:h 


h(s,fi,e) 


+  AL, 
1 
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By  equation  (2.2),  we  know  that 
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where  K3(tf )  =  tf  rriax{  T\>tj  exp(Kitf),  1}.  By  Gronwall’s  inequality,  we  then  have  that 

h(t,n,e)\  rl 


<  K3(tf)  {|Ai,jf|  +  |A^f|  +  |Ay  }d8'eK^. 


>(t)  ~ 

Since  |A4/jr  |,  |Ay  |,  |  A|,  s|  — >  0  uniformly  in  s'  as  |e|  — >  0,  we  can  then  conclude  that  for  t  £  [— r,  tf],  h/e  — >  v  as 
|c  |  — »■  0.  Therefore,  the  partial  derivative  of  x  with  respect  to  ji  i  (evaluated  at  =  g  >  0)  exists  and  satisfies  (2.1)  with 
the  initial  condition  (^(0),  iP)  =  (0,  0)  G  l4  x  cif(—r,  0;  M4),  which  completes  the  proof. 


Remark  2.4.  The  line  of  reasoning  presented  here  in  Lemma  2.1  and  Theorem  2.3  concerns  the  existence  and  continuity 
(in  t)  of  the  derivative  of  a  solution  to  (1.1)  with  respect  to  the  specific  parameter  // 1 .  Similar  arguments  (with  minimal 
changes  to  g3)  establish  the  existence  and  continuity  (in  t)  of  derivatives  with  respect  to  g  o,  <ti  ,  and  a 2-  For  the  parameters 
that  appear  in  (1.1)  as  linear  coefficients,  g  i  and  r/j  are  slightly  altered  (dependent  upon  the  parameter  under  consideration), 
while  (73  =  0.  However,  these  differences  do  not  change  the  conclusion  that  the  derivative  of  the  solution  x{t)  (with  respect 
to  any  parameter  appearing  on  the  right  side  of  (1.1))  exists  and  is  continuous  in  time.  One  can  also  establish  differentiability 
of  solutions  with  respect  to  discrete  delays  (i.e.,  when  P\  or  A  is  a  Dirac  measure)  and  well-posedness  of  the  appropriate 
sensitivity  equations.  The  arguments,  while  in  the  spirit  of  those  given  above,  are  however  somewhat  more  tedious  and  will 
not  be  given  here. 
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3.  Analysis  and  Results 


In  this  section  we  examine  some  applications  of  the  theory  developed  in  Section  2.  All  of  the  simulations  presented  in 
this  section  were  done  using  Matlab  software  originally  developed  in  [3]  for  the  purpose  of  simulating  systems  of  Abstract 
Evolution  Equations  (AEE’s)  that  are  linear  in  the  delay  (e.g.,  system  (1.1)).  As  can  be  inferred  from  equation  (2.1),  in  order 
to  solve  sensitivity  equations,  one  needs  the  solution  x  of  the  original  system.  Therefore,  we  use  the  calculated  solution 
from  [3]  with  parameters  q*  G  (}n,i  (the  space  of  admissible  parameters)  identified  by  minimizing  the  cost 


(3.1) 


J(q)  = 


10 


\ 


10 


'£i(X(ti-,q)-Xi)\ 


i=  l 


over  q  £  Qaj  and  where  X  =  {A^,}  is  the  data  from  [40]  taken  at  times  {i,}.  However,  we  are  not  able  to  compute  the 
exact  solution  x  and  thus  (as  described  in  [3,  6])  we  minimize 


(3.2) 


J  {q)  =  w 
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£( XN(U;  v)  -  v,)" 


i= 1 


where  x  A  =  (VN ,  AN ,  CA  ,  TA )  is  an  appropriate  approximation  to  x  and  N  is  an  integer  describing  the  accuracy  of  the 
numerical  simulation  such  that  lim n^oc  xn  ( t ;  q)  =  x (t;  q).  The  numerical  scheme  (also  described  in  [3,  6])  is  such  that  as 
N  — >  oo,  a  minimizer  <jA*  to  (8)  converges  to  q* ,  a  minimizerto  (8).  Both  the  experimental  results  and  the  numerical  best 
fit  solution  xN  (using  parameters  gA  *  from  Table  2)  are  depicted  in  Figure  3.1. 


Parameter 

Value 

Units 

n. 4 

0.112 

hours-1 

TiC 

0.011 

hours-1 

7 

9E  -4 

hours-1 

<5.4 

0.078 

hours-1 

Sc 

0.025 

hours-1 

K. 

0.017 

hours-1 

8 

IE  -  12 

(cell  hours)-1 

P 

1.3E-6 

(cell  hours)-1 

Pi 

-22.8 

hours 

P-2 

-26 

hours 

Table  2.  Optimal  in  vitro  model  parameter  values. 


By  Theorem  2.3,  we  can  legitimately  consider  the  derivative  of  both  sides  of  (1.1)  with  respect  to  any  appropriate 
parameter.  We  first  consider  the  derivative  of  x(t)  with  respect  to  /i  i  at  /i  \  =  p 

+  ikhit)  forO  <t  <tf 


(3.3) 


■j^(x(0,/j,),xo  (fJ,)) 


^($(0),$)  £K4  x^(-r,0;K4). 
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If  we  denote  v(t)  =  //)  (for  some  specific  value  of  pi  =  p  >  0),  we  obtain  the  sensitivity  equations 

v(t)  =  gi(x(t]fj,)-,v(t))+g2(ii\Vt}+g3(xt(iJ,),ir,l)  forO  <t<tf 

(3.4) 

(i>(0),i>o)  =  (0,0)  £l4  x^(-r,0;l4), 

where  gi,  go,  g3  are  as  defined  in  Lemma  2.1.  As  before,  due  to  the  complexity  of  the  right  side  of  (3.4),  we  cannot 
solve  exactly  for  the  solution  v(t).  Moreover,  we  do  not  have  x  which  appears  in  the  terms  g  t  and  g3 ;  we  only  have 
an  approximation  xN  to  x.  Therefore,  we  must  propose  a  viable  numerical  scheme  to  calculate  an  approximation  v  A  to 
solutions  of  (3.4),  with  x  replaced  by  xN  such  that  limjv-^op  vN  =  v. 

Hence  we  consider ; > A  an  approximate  solution  to  (3.4)  with  x  =  x  A  in  the  coefficients.  This  is  a  linear  nonautonomous 
system  of  the  form 

vN ( t )  =  £/N(t)vN ( t )  +  g-2{vt  )  +  gs{xt)  for  0  <  t  <  tf 

(3.5) 

(t^(0),^)  =  (0,0)  G  M4  x  Sf(-r,  0;M4) , 

where  N  is  fixed,  xN  is  given,  and  ,c/A  maps  M4  x  ^(—r,  0;  M4)  to  M4.  Note  that  this  is  a  special  case  of  the  systems 
treated  in  [3],  where  existence  and  uniqueness  are  guaranteed.  To  obtain  convergence  of  v  A  to  v  (the  unique  solution  to 
(3.4)),  we  turn  to  [2],  A  straightforward  extension  of  the  theory  presented  in  [2]  to  treat  nonautonomous  linear  systems  such 
as  (3.5)  will  yield,  (under  the  approximation  scheme  described  in  [3]),  the  desired  convergence. 

If  we  were  to  plot  simulations  of  (3.4)  (or  actually,  the  approximate  solutions  defined  by  (3.5)),  interpretations  of  these 
plots  would  suggest  specific  effects  that  changes  in  // 1  would  have  on  the  solution  x.  Moreover,  if  we  were  to  also  perform 
the  analogous  derivation  for  the  infection  rate  p,  a  plot  of  that  sensitivity  function  would  depict  the  effect  that  changes  in  p 


Figure  3.1.  Data  from  [40]  and  best  fit  simulation  xN  of  (1.1). 
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would  have  on  x.  Since  /i  \  and  p  differ  in  their  units,  the  sensitivity  functions  for  p  \  and  p  would  also  have  different  units, 
thus  rendering  any  comparison  meaningless.  We  turn  to  the  sensitivity  analysis  literature  to  resolve  this  issue.  To  enable 
a  comparison  of  the  effects  that  parameters  with  different  units  have  on  the  solution,  we  simply  multiply  by  the  parameter 
under  consideration,  e.g.,  (  (f;  p),  t^;eo  (f;  p),  (i;  yu) ,  (f;  /r.))  •  p.  This  form  of  the  sensitivity  function  is 

known  as  the  semirelative  or  semilogarithmic  or  unnormalized  sensitivity  function  [10,  1 1].  Moreover,  this  form  is  actually 
the  differential  of  x  with  respect  to  p  i  at  p  in  the  direction  p 

d 

D^Xi{t-,n)[n]  =  (- — Xj(t;  n))  ■  p 
a 

for  i  =  1,  2, 3, 4.  With  this  weighting,  we  now  have  the  tools  to  rank  the  parameters  with  regard  to  their  influence  over  the 
solution. 

Figure  3.2  depicts  the  approximation  vN  of  the  solution  v  to  the  (3.4)  (at  p  =  —22.8),  with  each  compartment  multiplied 
by  /i.  It  is  important  to  realize  that  while  the  y-axis  in  Figure  3.2  has  units  of  cells  or  virions  respectively,  it  should  still 
be  thought  of  as  a  plot  reflecting  changes  in  the  state  with  respect  to  changes  in  p  l ■  I'1  other  words,  we  interpret  the 
upper-left  plot  of  Figure  3.2  to  suggest  that  for  a  (positive)  change  in  the  mean  delay,  the  virion  compartment  V  will  be 
dramatically  smaller  just  before  day  10  and  then  larger  around  day  12  (relative  to  V(t;  —22.8)).  Likewise  for  a  change  in 
Hi,  the  acutely  infected  cell  compartment  A  will  be  slightly  smaller  around  day  9  and  dramatically  larger  around  day  10 
(relative  to  A(t ;  —22.8)).  All  the  plots  depicted  in  Figure  3.2  suggest  that  there  will  be  dramatic  changes  in  the  solution  for 
changes  in  p\ ,  and  indeed  Figure  3.3  supports  this  claim  (as  well  as  the  specific  predictions  suggested  by  the  interpretation 
of  Figure  3.2).  For  this  simulation,  it  is  important  to  note  that  there  is  practically  no  indication  that  the  solution  x  will  exhibit 
any  sensitive  to  hi  until  around  day  5.  In  other  words,  for  simulations  on  a  short  time  interval  (i.e.,  t  £  [—r,  120]  hours), 
one  could  easily  conclude  that  the  solution  x  is  insensitive  to  p  i  (in  the  neighborhood  of  hi  =  H  =  —22.8  hours). 

As  another  example,  let  us  consider  the  solution  parameterized  with  respect  to  the  infection  rate  p,  i.e.,  x(t)  =  x(t;p). 
Thus  the  derivative  of  ( 1 . 1 )  with  respect  to  p  at  p  =  1.3  x  1 0  " 6  is 


(3.6) 


J^,x{t\p)  =  -^L{x{t]p%xt{p))  +  -^f1{x{t]p),p)  +  -^f2{t)  forO  <t<tf 

£(x(0,p),xo(p))  =  £(*(0),*)  GK4  x^(-r,0;IR4). 


As  mentioned  in  the  last  part  of  Section  2,  the  sensitivity  equations  with  respect  to  different  parameters  will  be  slightly 
different  than  (3.4),  but  unique  solutions  still  exist  and  are  continuous  (for  each  system  of  sensitivity  equations).  Figures  3.5 
and  3.4  depict  the  semirelative  sensitivity  functions  for  p  and  p  -j,  respectively  .  A  comparison  of  the  scales  on  the  vertical 
axis  in  Figure  3.2  versus  the  axis  in  Figures  3.5  and  3.4  suggests  that  changes  in  p  l  have  a  more  significant  influence  in  the 
solution  x  than  changes  in  p  >  or  p  (and  in  one  of  the  compartments  by  over  four  orders  of  magnitude).  This  result  coincides 
nicely  with  one  of  the  primary  conclusions  from  [3],  in  which  we  concluded  that  when  fitting  the  data,  adding  the  second 
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Days 


Figure  3.2.  Simulation  of  the  semirelative  sensitivity  solution  with  respect  to  /r,  i  at  ji  \  =  /. 


x  10 


Virions  V 


x  ^  q6  Acutely  Infected  Cells  A 


x  o6  Uninfected  Target  Cells  T 


Figure  3.3.  Simulations  of  (t;  —24.8)  and  xN  (t;  —22.8). 
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delay  between  than  acute  and  chronic  infection  was  not  as  significant  as  inclusion  of  the  delay  between  viral  infection  and 
viral  production. 


0  5  10  15  20  25 

Days 


x  IQ5  Chronically  Infected  Cells  C 


Days 


x  io6  Uninfected  Target  Cells  T 


Figure  3.4.  Simulation  of  semirelative  sensitivity  solution  with  respect  to  the  infection  rate  p  for  p  = 
1.3  x  106. 


Now  that  we  have  established  the  framework  for  calculating  semirelative  sensitivity  functions,  let  us  consider  how  to 
rank  the  influence  that  changes  in  the  individual  parameters  have  upon  the  solution  x.  Clearly,  there  are  many  options,  but 
for  simplicity,  we  will  rank  the  parameters  according  to  the  magnitude  of  the  oo-norm,  e.g.,  for  the  virion  compartment  and 
the  sensitivity  with  respect  to  6  a,  we  consider 

max  \Ds,V{t-,  0.0776)  [0.077611  . 
te[o,*/] 

To  illustrate  our  reasoning,  we  will  focus  on  just  the  virion  compartment  V .  Of  the  parameters  over  which  we  performed  our 
NLS  in  [3],  the  chosen  metric  was  largest  for  the  parameters  pi,  ha,  6a,  and  6U.  Figure  3.6  depicts  (for  the  compartment 
V’),  the  absolute  values  of  the  semirelative  sensitivity  functions  with  respect  to  p  i ,  ha,  6a,  and  6U,  for  t  G  [8.5, 15]  (the 
domain  where  there  is  the  most  activity  in  the  sensitivity  functions).  The  interpretation  of  this  figure  strongly  suggests  that 
6a  and  n  \  have  the  strongest  influence  over  the  solution  in  the  virion  compartment  (in  the  chosen  oo-norm).  Therefore,  for 
the  use  of  equation  (1.1)  (as  a  model  to  simulate  HIV  pathogenesis),  both  the  viral  production  rate  and  the  death  rate  for 
acutely  infected  cells  (ua  and  6a  respectively)  should  be  given  top  priority  when  choosing  which  parameters  to  determine 
with  a  high  degree  of  accuracy.  In  other  words,  these  parameters  play  an  important  role  in  the  model  and  obtaining  good 
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values  for  them  is  more  important  to  the  system  response  than  good  values  for  other  parameters  to  which  solutions  are  less 
sensitive. 


Virions  V 


Days 


x  i  o4  Acutely  Infected  Cells  A 


Days 


Figure  3.5.  Simulation  of  semirelative  sensitivity  solution  with  respect  to  the  mean  delay  between  acute 
and  chronic  infection  /ij  for  fi  >  =  —26. 


Figure  3.6.  Absolute  value  of  simulations  of  semirelative  sensitivity  solutions  for  several  parameters 
(V  compartment  only). 
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4.  Concluding  Remarks 

As  discussed  in  Section  1.1,  the  taking  of  a  derivative  (with  respect  to  parameters)  of  the  equations  governing  a  system 
is  not  a  new  idea  and  indeed  has  been  around  (in  some  form)  for  at  least  170  years.  Within  control  theory  and  engineering 
applied  to  physical  systems,  the  forms  of  the  fundamental  mathematical  models  often  are,  for  the  most  part,  relatively  well 
established  and  not  so  open  to  debate.  For  example,  in  some  investigations,  it  may  not  be  fruitful  to  question  the  significance 
of  the  viscosity  parameter  in  the  Navier-Stokes  equations  (although  sensitivity  of  flow  patterns  to  viscosity  is  sometimes  very 
important,  see  [42]).  However,  the  constitutive  parameters  and  forms  of  the  mathematical  models  employed  in  the  biological 
sciences  are  frequently  not  as  well  agreed  upon,  and  indeed  (as  is  evidenced  by  the  literature)  open  to  considerable  debate. 
Since  the  current  approach  to  sensitivity  was  originally  developed  in  the  context  of  control  theory,  the  cited  literature  is 
(understandably)  biased  toward  that  field;  a  considerable  proportion  of  the  papers  are  devoted  to  analyzing  the  sensitivity 
of  transfer  functions  and  eigenvalues.  Thus  the  application  of  mathematically  rigorous  sensitivity  analyses  to  dynamical 
systems  designed  to  model  biological  phenomena  does  not  seem  to  be  common  practice.  Indeed,  many  sensitivity  studies 
often  involve  copious  simulations.  As  such,  there  are  many  possibilities  that  have  not  been  fully  examined. 

In  the  analysis  presented  in  this  paper,  we  only  considered  first  derivatives  of  the  components.  In  theory,  we  could 

o  2 

have  examined  derivatives  with  respect  to  multiple  parameters  (joint  sensitivities),  an  analysis  of  which  could  be 

used  to  ascertain  the  independent  identifiability  of  parameters.  We  could  have  also  taken  a  derivative  with  respect  to  the 
initial  conditions,  which  (as  is  intuitive)  would  suggest  the  influence  of  the  initial  conditions  over  the  solution  (this  can  be 
extremely  useful  in  certain  biological  investigations).  Finally,  we  could  have  considered  the  derivative  of  the  least  squares 
functional  (3.1)with  respect  to  a  parameter  (as  was  explored  in  [22]),  which  could  then  be  used  as  part  of  a  jacobian  in  an 
optimization  algorithm  (as  part  of  a  parameter  estimation  scheme). 

The  process  of  taking  the  derivative  of  a  system  with  respect  to  a  parameter  is  usually  not  an  exceedingly  challenging  task 
and  it  is  important  to  remember  that  the  sensitivity  function  only  reveals  the  local  behavior  (since  it  is  a  derivative)  around 
the  fixed  parameter  value.  However,  this  idea  can  yield  useful  insights  into  the  solution  of  complex  systems  (even  those 
with  nonlinearities  and  delays)  such  as  (1.1).  Effectively,  the  technique  of  using  simulation  sensitivity  functions  presented 
in  this  paper  is  a  more  efficient  (and  mathematically  rigorous)  way  to  attain  insight  into  a  system  than  manually  adjusting  a 
parameter  and  observing  the  effect  on  the  solution  through  massive  simulation  efforts. 
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